knitr::opts_chunk$set(warning = FALSE, message = FALSE)
options(warn = -1)
setwd("C:/Users/Angelica/Documents/Movilidad/Data")
library(readxl)
Hogares <- read_excel("a. Modulo hogares.xlsx")
#Vehiculos <- read_excel("b. Modulo vehiculos.xlsx")
#Personas <- read_excel("c. Modulo personas.xlsx")
#Viajes <- read_excel("d. Modulo viajes.xlsx")
#Etapas <- read_excel("e. Modulo etapas.xlsx")

ESTADÍSTICAS DESCRIPTIVAS DE HOGARES

UTAM: Unidades Territoriales de Análisis de Movilidad

Cargamos el archivo shapefile

library(sf)
library(ggplot2)
setwd("C:/Users/Angelica/Documents/Movilidad/UTAM")
UTAM_shape <- st_read("UTAM2023.shp")
## Reading layer `UTAM2023' from data source 
##   `C:\Users\Angelica\Documents\Movilidad\UTAM\UTAM2023.shp' using driver `ESRI Shapefile'
## Simple feature collection with 142 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.44978 ymin: 3.73103 xmax: -73.86269 ymax: 5.041663
## Geodetic CRS:  MAGNA-SIRGAS

Calculamos el conteo de hogares

El codigo UTAM del shapefile es diferente al de hogares dado que, por ejemplo, en hogares la UTAM1 se escribe UTAM001. Por tanto es necesario modificar la variable “cod_utam_hg”

library(dplyr)
Hogares <- Hogares %>%
  mutate(cod_utam = gsub("^(UTAM|UPR)0*", "\\1", cod_utam_hg))

num_hogares <- as.data.frame(table(Hogares$cod_utam))
colnames(num_hogares) <- c("UTAM", "Hogares")  

Luego, lo unimos con el shapefile y generamos el mapa

shape_joined <- UTAM_shape %>%
  left_join(num_hogares, by = "UTAM")

ggplot() +
  geom_sf(data = shape_joined, aes(fill = Hogares), color = "black") +
  scale_fill_viridis_c(option = "C", na.value = "gray90") +
  theme_minimal() +
  labs(title = "Número de hogares por UTAM", fill = "Hogares")

Si se excluyen las zonas de planeamiento rural (UPR) podemos ver más detalladamente las UTAM centrales

shapefile_utam <- shape_joined %>%
  filter(grepl("^UTAM", UTAM))

ggplot() +
  geom_sf(data = shapefile_utam, aes(fill = Hogares), color = "black")+
  scale_fill_viridis_c(option = "plasma", na.value = "grey90") +
  theme_minimal() +
  labs(title = "Número de hogares por UTAM",
       fill = "Hogares")

Para solucionar las dimensiones del mapa, se puede hacer un gráfico interactivo con plotly

library(stringr)
shapefile_utam <- shape_joined %>%
 mutate(UTAMNombre = if_else(str_trim(str_to_upper(UTAMNombre)) == "N/A", MUNNombre, UTAMNombre))%>%
 mutate(label = paste0("Codigo : ", UTAM, "\nNombre UTAM: ", UTAMNombre, "\nHogares: ", Hogares ))
  
mapa <- ggplot() +
  geom_sf(data = shapefile_utam, aes(fill = Hogares, text = label), color = "black") +
  scale_fill_viridis_c(option = "plasma", na.value = "grey90") +
  theme_minimal() +
  labs(title = "Número de hogares por UTAM", fill = "Hogares")

library(plotly)
ggplotly(mapa, tooltip = "text")

Visto en diagrama de barras para municipios diferentes a Bogotá

library(forcats)
ggplot(subset(Hogares, nom_mpio_hg != "Bogotá D.C."), 
       aes(x = fct_rev(fct_infreq(nom_mpio_hg)))) +
  geom_text(stat = "count", aes(label = ..count..), hjust = -0.2, size = 3) +
  geom_bar(fill = "steelblue") +
  coord_flip() +
  labs(
    title = "Número de hogares por municipio",
    x = "Municipio",
    y = "Cantidad de hogares"
  ) +
  theme_bw()

Variable ingresos

ggplot(Hogares, aes(x = ingre_mes_hg)) +
  geom_bar(fill = "steelblue") +
  geom_text(stat = "count", aes(label = ..count..), vjust = -0.3, size = 3) +
  labs(
    title = "Distribución de ingreso mensual",
    x = "Ingreso mensual (COP)",
    y = "Número de hogares"
  ) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

# Contar hogares por municipio e ingreso
tabla_prop <- Hogares %>%
  filter(ingre_mes_hg != "NS/NR") %>%
  count(nom_mpio_hg, ingre_mes_hg) %>%
  group_by(nom_mpio_hg) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup()

# Heatmap con proporciones
ggplot(tabla_prop, aes(x = ingre_mes_hg, y = fct_infreq(nom_mpio_hg), fill = prop)) +
  geom_tile(color = "white") +
  scale_fill_viridis_c(option = "C", labels = scales::percent) +
  labs(
    title = "Distribución de ingreso mensual por municipio",
    x = "Ingreso mensual",
    y = "Municipio",
    fill = "Porcentaje"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

ggplot(subset(Hogares, nom_mpio_hg == "Bogotá D.C."),
       aes(x = estrato_hg, fill = factor(ingre_mes_hg))) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = colores) +  # Usamos scale_fill_manual con los colores generados
  labs(
    title = "Distribución de ingreso mensual por estrato (Bogotá)",
    x = "Estrato",
    y = "",
    fill = "Ingreso mensual"
  ) +
  theme_bw()

Quitando a categoría NS/NR

  hogares_bog <- Hogares %>%
  filter(nom_mpio_hg == "Bogotá D.C.", ingre_mes_hg != "NS/NR")

hogares_bog %>%  
  ggplot(
    aes(x = estrato_hg, fill = factor(ingre_mes_hg))) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = colores) +  # Usamos scale_fill_manual con los colores generados
  labs(
    title = "Distribución de ingreso mensual por estrato (Bogotá)",
    x = "Estrato",
    y = "",
    fill = "Ingreso mensual"
  ) +
  theme_bw()

O visto no en proporcion

ggplot(subset(hogares_bog, nom_mpio_hg == "Bogotá D.C."),
       aes(x = ingre_mes_hg, fill = factor(estrato_hg))) +
  geom_bar(position = "dodge") +  # Posición dodge para que las barras estén separadas
  labs(
    title = "Distribución de ingreso mensual por estrato (Bogotá)",
    x = "Estrato",
    y = "Número de hogares",
    fill = "Ingreso mensual"
  ) +
  theme_bw()

Tamaño de los hogares

ggplot(Hogares, aes(x = perstotal_hg)) +
  geom_histogram(binwidth = 1, fill = "steelblue", color = "white") +
  labs(
    title = "Distribución tamaño del hogar",
    x = "Número de personas en el hogar",
    y = ""
  ) +
  theme_bw()

ESTADÍSTICAS DESCRIPTIVAS DE VEHICULOS

Tipos de vehículo

Se cuenta con una categorización bastante amplia con frecuencias bajas para algunos tipos de vehículo

kable(tabla_vehiculos, caption = "Tipos de vehículo")
Tipos de vehículo
Tipo de Vehículo Cantidad
Automóvil 6441
Automóvil o camioneta de servicio especial 64
Bicicleta motor eléctrico 55
Bicicleta motor gasolina 9
Bicicleta niños 329
Bicicleta sin motor 4691
Bicitaxi motor eléctrico 9
Bicitaxi motor gasolina 13
Bicitaxi sin motor 190
Camión 83
Campero/Camioneta 901
Moto - carro 134
Motocicleta 4028
Otro 69
Patineta 70
Patineta eléctrica 92
Pick Up/Van 69
Taxi 115
Triciclo - moto 3
Vehículo tracción animal 3
Vehículo tracción humana 24

Por tanto, es necesario realizar una recategorización de la variable. Se propone la siguiente manera

Recategorización de tipos de vehículo
Categoría.nueva Tipos.originales.incluidos
Automóvil Automóvil
Motocicleta Motocicleta
Campero/Camioneta Campero/Camioneta
Otros motorizados Taxi, Automóvil o camioneta de servicio especial, Moto - carro, Camión, Pick Up/Van, Triciclo - moto, Bicicleta motor gasolina, Bicitaxi motor gasolina, Bicicleta motor eléctrico, Bicitaxi motor eléctrico, Patineta eléctrica
Otros Bicicleta sin motor, Bicicleta niños, Bicitaxi sin motor, Patineta, Vehículo tracción humana, Vehículo tracción animal
ggplot(Vehiculos, aes(x = fct_rev(fct_infreq(tipo)))) +
  geom_bar(aes(fill = tipo)) +
  labs(
    title = "Distribución tipos de vehículos",
    x = "",
    y = "Cantidad de vehiculos"
  ) +
  theme_bw() +
  coord_flip() +
  scale_fill_brewer(palette = "Set2") +
  theme(legend.position = "none")

Tipo de combustible

Se propone una recategorización de la variable

Recategorización tipo de combustible
Categoría nueva Categoria Original
Gasolina Sólo Gasolina
Diésel Diésel
Eléctrico Eléctrico
No aplica No aplica
Otro GNV, GNV y gasolina, Híbrido (eléctrico-gasolina / diésel)
ggplot(Vehiculos, aes(x = fct_rev(fct_infreq(combustible)))) +
  geom_bar(aes(fill = combustible)) +
  geom_text(stat = "count", aes(label = ..count..), hjust = -0.2, size = 3) +
  labs(title = "Distribución por tipo de combustible",
       x = "Tipo de combustible", y = "Cantidad de vehiculos")+
  theme_bw() +
  coord_flip() +
  scale_fill_brewer(palette = "Set2") +
  theme(legend.position = "none")

ggplot(Vehiculos, aes(x = (fct_infreq(tipo_placa)))) +
  geom_bar(aes(fill = tipo_placa)) +
  labs(title = "Distribución por tipo de placa",
       x = "Tipo de placa", y = "Cantidad de vehiculos")+ 
theme_bw() +
  scale_fill_brewer(palette = "Set2") +
  theme(legend.position = "none")

UTAM

Grafico excluyendo zonas de planeamiento rural UPR

Vehiculos %>%
  ggplot(aes(x = factor(estrato_hg))) +
  geom_bar(aes(fill = tipo),position = "fill") +  
  scale_y_continuous(labels = percent_format()) +
  labs(
    title = "Distribución del tipo de vehículo por estrato",
    x = "Estrato",
    y = "Proporción",
    fill = "Tipo de vehículo"
  ) +
theme_bw()

library(tidyr)
#Filtrar Vehículos motorizados
vehiculos_filtrados <- Vehiculos %>%
  filter(tipo != "Otros")

#Contar vehículos por hogar
hog_vehiculos <- vehiculos_filtrados %>%
  group_by(cod_hog) %>%
  summarise(num_vehiculos = n())

#Unir con Hogares y reemplazar NA por 0
Hogares_vehiculos <- Hogares %>%
  left_join(hog_vehiculos, by = "cod_hog") %>%
  mutate(num_vehiculos = replace_na(num_vehiculos, 0))

ggplot(Hogares_vehiculos, aes(x = factor(num_vehiculos))) +
  geom_bar(fill = "steelblue") +
  labs(
    title = "Distribución de número de vehículos motorizados en el hogar",
    x = "Cantidad de vehículos",
    y = "Número de hogares"
  ) +
  theme_bw()

Creo que en hogares cuyo ingreso es de 0-0.4M debe haber bicitaxis que generan ese incremento en el numero de vehiculos porque no me parece normal